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The square- gradient density-functional model with triple-parabolic free energy, that was used 
previously to study the homogeneous bubble nucleation [J. Chem. Phys. 129, 104508 (2008)], is 
used to study the stability of the critical bubble nucleated within the bulk under-saturated stretched 
fluid. The stability of the bubble is studied by solving the Schrodinger equation for the fluctuation. 
The negative eigenvalue corresponds to the unstable growing mode of the fluctuation. Our results 
show that there is only one negative eigenvalue whose eigenfunction represents the fluctuation that 
corresponds to the isotropically growing or shrinking nucleus. In particular, this negative eigenvalue 
survives up to the spinodal point. Therefore the critical bubble is not fractal or ramified near the 
spinodal. 

PACS numbers: 47.55.db, 64.60.qe, 82.60.Nh 
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I. INTRODUCTION 

The stability of the liquid- vapor interface and its rela- 
tion to correlations and capillary waves has been studied 
for more than four decades^"— . The stability of the in- 
terface is determined from the eigenvalue problem of the 
stability matrix—iS or the Schrodinger equatior^'^'^'i de- 
rived from the classical density functional theory^,. It is 
well recognized that there is always zero eigenvalue which 
corresponds to the free translation of the planar liquid- 
vapor interface^. The detailed study of the eigenvalues 
of the stability matrix and the Schrodeinger equation has 
already been made by Bukman et al^ for the planar inter- 
face and by Varea and Roblcdc^ for the curved interface. 

Very recently, a renewed interest in the stability of the 
liquid-vapor interface has been revived^"—. Uline and 
Corti^'^° have studied the stability of the spherical liquid- 
vapor interface of the critical bubbleSs- and droplet^° of 
the liquid-vapor nucleation, and has shed doubt on the 
classical view of the so-called classical nucleation the- 
orjiiiii^. Their work has induced debates^"— on the 
validity of the concept of the minimum-free energy path 
on the free-energy landscape of the nucleation. 

In this paper, we will use a simple square-gradient 
density functional theory with a triple-parabolic free en- 
ergy proposed by Granasy and Oxtobj^, which has been 
used to study the various properties of the critical bub- 
ble of homogeneous bubble nucleation by the author—, 
to study the stability of the liquid- vapor interface of the 
critical bubble. We choose this square-gradient density 
functional model as it captures the most basic properties 
of nucleation, yet many physical quantities can be han- 
dled analytically. In fact, very recently, Li and Wilem- 
ski^S have compared the results obtained from the accu- 
rate density functional theory with the result from the 
approximate square-gradient theory and found that the 
two results agree qualitatively well. In particular, the 



eigenvalue problem of the stability matrix reduces to the 
solution of the standard Schrodinger equation for a parti- 
cle in a potential well^'^'^'i in the square-gradient density 
functional model which is much simpler than the eigen- 
value problem of the stability matris^i^iiS that consists 
of roughly 10^ to 10^ elements in the density functional 
model. 

Although Varea and Robledoi used a similar square- 
gradient model with a triple-parabolic free energy to 
study the stability of curved interface, they were inter- 
ested in the stability of curved interface in general and 
paied less attention to the nucleation problem. In this 
paper, we will use a more realistic triple-parabolic free 
energyiS. and pay most attention to the stability of the 
critical nucleus of homogeneous nucleation including the 
nucleus near the spinodal point. 

In Section II of this paper we will present a short review 
of the stability of the spherical liquid- vapor interface in 
the square-gradient density-functional model. In Section 
III, we will present the numerical results for the stability 
of the spherical liquid- vapor interface of the critical bub- 
ble within the triple-parabolic free energyi^ and discuss 
the implications of the results in light of the stability of 
the spherical critical nucleus against non-spherical fluc- 
tuations. Finally Section IV is devoted to the concluding 
summation. 



II. STABILITY OF THE LIQUID- VAPOR 
INTERFACE IN THE SQUARE-GRADIENT 
DENSITY-FUNCTIONAL MODEL 



In the square-gradient density-functional model of the 
fluid^ '^^i^^'^^ , the free energy (grand potential) W of the 
inhomogeneous fluid, such as the critical bubble in the 
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stretched liquid is written as 



(1) 



where c is the coefficient of square gradient, and Aw ((/)) 
is the local grand potential density as the function of the 
local density <^(r). This form of the free energy is also 
known as the Cahn-Hilliard model^^ or the phase-field 
model^i^. 

The stationary density profile (ps of the critical nucleus 
can be determined from the stationary condition of the 
grand potential written as 



5W [<t>{v)]/5<t>{v)^Q, 
which leads to the Euler-Lagrange equation 



d(l) 



2cV^0 (r) = 0. 



(2) 



(3) 



Now the stationary profile 0s (r) and the work of forma- 
tion Ws of the critical nucleus can be obtained by solving 
the differential equation Eq. 

The stability of this stationary profile 0s will be stud- 
ied from the second variation of the functional W ^ which, 
after integration by parts, becomes 




5(l> (r) - 2cV2,50 (r) ) (fr, 
(4) 

where Scf) (r) is a small variation of the order parameter 
from the stationary profile 0s (r). Let ipn be the eigen 
functions of the Schrodinger equatior^ii^i^ 



2cVV« (r) + V (r) 7/;„ (r) = S„V« (r) , (5) 



where 



v{r) 



plays the role of the potential and E„ is the eigenvalue. 
Then the order-parameter fluctuation can be expanded 
as 



50 (r) = ^C„V„ (r), 

n 

and the second variation takes the form 



(7) 



(8) 



The stability of the stationary solution 0s (r) will be 
studied from the sign of the eigenvalues E,-i}-7 . In fact, by 
using the non-conserved order-parameter dynamics^^'^ 



90 _ 5W 
'dt ~^ 64) (r) ' 



(9) 



the time evolution of the order parameter (r, t) near 
the stationary critical nucleus 0s is written as 

0(r,i)-0s(r)+<50(r,i), (10) 

and the time evolution of the fluctuation is given by^ 



. (r, t)^Yl (^) «^^P i-^Ent) , 



(11) 



where F characterizes the time scale of the evolution. 
The negative eigenvalue En < indicates the unsta- 
ble growing mode of the order-parameter fluctuation^. 
Therefore, the bound state of the Schrodigner equation 
^ with negative eigenvalue i?„ < corresponds to the 
unbound growing mode of the fluctuation. Apparently, 
the stability analysis is meaningful only near the station- 
ary critical nucleus which satisfies Eq. 

The stability of critical nucleus can be studied by solv- 
ing the simple one-body problem of the Schrodinger equa- 
tion ([U in the simple square-gradient density-functional 
model. In contrast, in the original density functional 
model®'^, the stability problem leads to the integral equa- 
tion, which is usually transformed into the eigenvalue 
problem of the stability matrix with roughly 10'^ to 10^ 
elements. 



III. STABILITY OF THE CRITICAL BUBBLE 

A. Triple-Parabolic Model Free Energy 

In order to obtain a more detailed description of 
the stability of the critical bubble, we chose the triple- 
parabolic model for the free energy Aw (0) in Eq. |T]) used 
previously to study the scaling properties of the critical 
bubble^^ originally proposed by Granasy and Oxtoby^*: 



(6) Aw(0) 



Ao 
2 

2 



Am, 



< 



A/i + Wo, 



A2 

2 



0yl < < 0B, 

(12) 

with Ao, A2 > and Ai < 0, which consists of three 
parabolas centered at the vapor density 0o, and at the 
free energy barrier 0i , and at the liquid density 02 , which 
we call "vapor", "spinodal" and "liquid" part of the free 
energy. 

The parabolic curvatures Ao and A2 are related to the 
compressibility of vapor and liquid phases^^, and A/i is 
the free energy difference between the liquid and the va- 
por. Since we can write A/x in Eq. \12\ by A/i0 ~ A/i02, 
A/i is, in fact, the reduced chemical potential of the fiuid 
divided by the liquid density 02. The terminology "over- 
saturation" is used when A/i is positive and "under- 
saturation" when A/i is negative. Since a stretched liquid 
will be considered in this study, we will be mainly con- 
cerned with the under-saturation when A/i < 0. From 
the continuity of the free energy Aw(0), the boundaries 
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(j)A and (pB, as well as the barrier height loq and its loca- 
tion are given as the functions of Aq, Ai, A2 and the 
under-saturation A/i^^. 

The liquid spinodal is defined when the metastable liq- 
uid phase at (j)2 becomes unstable. This is realized when 
4'2 — 4>i, which leads to the under-saturation for the liq- 
uid spinodali^, 



spin 



1 AolAil 
'2Ao-f |Ai 



(02 - 00 ) 



(13) 



of the stretched liquid. 

In contrast to the previous models22r— where the com- 
pressibility diverges continuously as the spinodal is ap- 
proached, the compressibility remains finite until the 
spinodal point is reached in our triple-parabolic model 
as the curvature Aq and A2 is fixed. 
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FIG. 1: (a) The triple-parabolic free energy from CNT regime 
near the coexistence to the spinodal regime near the liquid 
spinodal for the case (i) (Table (b) The corresponding 
critical bubble at the CNT regime (A/i/A/ispin = 0.3) and at 
the spinodal regime (A/i/A/ispin = 0.7). 

In Fig. [TJa) we show the typical shapes of the triple- 
parabolic free energy Aw (0)^^. The corresponding den- 
sity profile 4>s {r) of the critical bubble is shown in 
Fig. [ijb). The radii ta and are the matching radius 
that satisfies 4>{rA) = and (t>{rB) = 4>b- Since the 
free energy consists of three parabolas corresponding to 
the vapor, spinodal and liquid parts, the density profile of 
the critical bubble consists of three parts that correspond 

0.3 



spin 



to the three parts of free energy when A^/A/i. 
near the coexistence. However, as the under-saturation 
A/z increases (|A^| becomes large) and it approaches the 



liquid spinodal A/igpin, the matching radius rA vanishes 
and the density profile consists of only two parts that cor- 
respond to the spinodal and the liquid parts (Fig. [TJb)). 
We use the terminology "CNT regime" for the former 
regime near the coexistence where the classical nucleation 
theory (CNT) is expected to be qualitatively correct, and 
" spinodal regime" for the latter near the spinodal where 
the spinodal nucleatio n^^'^^" — is expected to occur. 



B. Stationary Profile of the Critical Bubble 

1. CNT regime 

Density profile of a spherically symmetric critical bub- 
ble can be obtained from the Euler-Lagrange equation 
Eq. dSI: 

1 d f ^d(i)\ 1 9 Aw 



dr \ dr J 2c 
which leads to the ordinary differential equations 

d2$ 



dr2 



2 d$, o 

-— ^±r2$, = 0, * = 0,1,2, 
r dr 



(14) 



(15) 



for the three parabolas in Eq. (|12p , where F^ = yj\\i\/2c 
and $i(r) = 0(r) — 0,;, and -f sign is used for i — \ and 
— is used for i = 0,2 for ±. These differential equa- 
tions should be solved with appropriate boundary condi- 
tions^^. 

The solutions of this Euler-Lagrange equation in 
Eq. ([15]) for the critical bubble are given by 

^o{r) = $oA»-ACSch(rorA)sinh(For)/r, 

$i(r) = csc(Fi(r^ -^b)) 

X (-^isrs sin (Fi (r - r^)) 
+^iArA sin (Fi (r - rB)))/r, 

$2(?') = $23^5 exp (-F2r + F2rB) /r, 

for the three parts i = 0, 1, 2 respectively, where 



(16) 





Pa ~ 00 




Pa - 01 


$is = <; 


tiB - 01 


*2S = 


0B - 02 



(17) 



Finally, the matching radii rA and rs are determined 
from the simultaneous equations 



d^2 




d^i 


dr 




dr 








dr 




dr 



(18) 



where only $1 is a function of both rA and rs (Eq. (fT6|)). 
These simultaneous equations can be solved numerically 
using standard algorithms such as the Newton-Raphson 
method. Even though multiple roots of Eq. (IT5|) exist, the 
smallest and should be chosen as they corresponds 
to the nucleus with the lowest free-energy. 
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2. Spinodal regime 

In this case ta becomes zero. Therefore i = CNT 
part (see Fig. [T|) of the free energy density Aw (0) in 
Eq. and its sohition $o disappears. Then, the solu- 
tion for the Euler-Lagrange equation for $2 is the same 
as Eq. but the sohition for $1 now reads 

$i(r) = $iBrBCSC (Eirs) sin (Eir) /r, (19) 

for the critical bubble in the spinodal regime. 

In this case, the matching radius tb is simply deter- 
mined from the equation 



dr 



dr 



(20) 



which is explicitly written as 

TB (AaEi cot (Eirs) - |Ai|r2) - A2 + |Ai 



(21) 

that does not depend on the under-saturation A/i. 
Therefore, the matching radius is constant in the 
spinodal regime. Again this equation should be solved 
numerically and the smallest radius tb should be cho- 
sen. 

In Fig. [Ijb), we showed the typical density profiles in 
the CNT regime and in the spinodal regime. The critical 
bubble is larger in the CNT regime than in the spin- 
odal regime. However, the density difference between 
the inside and the outside of the bubble becomes smaller 
in the spinodal regime than in the CNT regime. Cor- 
respondingly, the interfacial thickness looks diffuse^^ as 
the spinodal is approached. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

A|j./A)l,p,„ 

FIG. 2; The matching radii rA and rs as functions of the 
scaled under-saturation A^/A/ispin for the case (i) and (iii) of 
the free energy parameter hsted in Table [T] In the spinodal 
regime, rA disappears, and r_B becomes constant. 

Figure [2] shows the matching radii and as func- 
tions of the scaled under-saturation A/i/A/Xgpin for the 
free energy parameter tabulated in Table U] The match- 
ing radius ta becomes zero as the under-saturation enters 
the spinodal regime, while the matching radius sur- 
vives and it becomes constant and independent of the 
under-saturation A/i in the spinodal regime as predicted 
from Eq. (pH . 



TABLE I: Three sets of the free energy parameters in Au {(p) 
of Eq. (|12p used in this work. Case (i) and (ii) were used 
previously^^ to check the universality of the scaling properties 
of the critical bubble. 



model 


c 


(po 


02 


Ao 


Ai 


A2 


Ao/IAil 


Ao/IAil 


case (i) 


1.0 


0.1 


1.0 


1.0 


-0.5 


0.3 


2.0 


0.6 


case (ii) 


1.0 


0.1 


1.0 


1.0 


-2.0 


0.9 


0.5 


0.45 


case (iii) 


1.0 


0.1 


1.0 


1.0 


-2.0 


1.4 


0.5 


0.7 



C. Eigenvalue problem 

The Schrodinger equation Eq. ([5]) for the triple- 
parabolic model becomes that for a particle in a three- 
dimensional square well potential: 



Ao, < r < rA, 
V (r) = - \Xi \ , rA <r <rB, 
A2, rB < r. 



(22) 



Figure 3 schematically shows the shape of the potential 
well. Depending on the magnitude of the potential bar- 
rier Ao and A2, there are two cases: (a) Aq > A2 and (b) 
A2 > Ao. 
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FIG. 3: The two cases of the square-well potential v (r) in 
Eq. (|22|l as the functions of the radial distance r. (a) Ao > A2, 
(b) A2 > Aq. The bound states have eigenvalues iS < A2. In 
case (b) when A2 > Ao, there always exist degenerate bound 
states Eni = Ao for I > 1. 

Since, we look for the negative eigenvalues, we are in- 
terested in the bound state whose eigenvalue E satisfies 



I All < E < X2 



(23) 



The eigenfunction ip (r) for the bound state of this spher- 
ically symmetric Schrodinger equation is written 



i'nlm (r) = Rnl (r) Yim {0 , (/)) , 



(24) 



using spherical polar coordinates {r,9,<j)), where 
Yim (9, (j)) is a spherical harmonic. The radial part of 
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the eigenfunction Rni (r) is determined from the radial 
equationiiii. 



j,2 ^j. 



,dR 
dr 



11 



R = 0. 



(25) 

Differentiating the Euler-Lagrange Eq. (fT4|) by r, one can 
easily show that there is always one bound state solution 
with E ^ and I = 1 whose eigenfunction is given by^ 

R (r) ^ d(t)/dr. (26) 

Other eigenvalues will be determined numerically. 

1. CNT regime 

In this regime, ta is finite and we have to solve 
the Schrodinger equation Eq. ([SS]) with three wells in 
Eq. When A2 < Aq (Fig. EJa)) the radial solution 

for the bound state with — |Ai| < i? < A2 is given by 

Aii{iaQr), r < va, 
Rni (r) = { ^ji ("1'') + C*"' ("1'^) ' rA <r <rB, 
Dhf'^ {ia2r) , tb <r, 



with A, B, C, and D are constant, and 



= v/(Ao - E) /2c, 
ai = ^{\Xi\+E) /2c, 
a2 = ^J[\2-E) /2c, 



(27) 



(28) 



where j;, and h^^'^ are the spherical Bessel, Neumann 
and Hankel functions, respectively^. 

The constants A, B, C and D are determined from 
continuity of the wave function Eq. (j27|) and its derivative 
that leads to the characteristic equation for the energy 
E, 



31 («0 -jiijl) -niirf) 
-J^ (C) -ni (C) 






/i« {^X) 
,(1)' 



(29) 



-O- (C) -C", (C) ^xh'i' i^x) 

where (x) = dji/dx etc. are the derivatives, and 

The roots of Eq. gives the eigenvalues En.ii 

0,1,2,...) for each angular momentum l{l = 0, 1, 2, 
Equation ([^ reduces, for example, to 



(30) 



1 



{e ^ [^cosh^ (-^cos {rj - () + xsin (77 - ()) 



- 77sinh$(xcos(77-C) + Csin(r7-C))]} = (31) 

for I — 0. Similar equations can be obtained for I > I. 
As has been noted in the previous subsection, we always 
have zero eigenvalues i?„=o,i=i = for / = 1. 



When A2 > Ao (Fig. [S^b)), Eqs. ^ to can be 
applicable as far as — |Ai| < i? < Aq. However, when the 
eigenvalue E falls within the range Xq < E < X2, the 
radial solution for r < in Eq. (|27p has to be replaced 
by 



Rni (r) = Aji (aor) , r < ta, 
and ao in Eq. (j28p is now given by 



ao = y/{E-Xo)/2c, 
and Eq. (|29|) becomes 




(1) 



(32) 



(33) 



ii'i iC) -vji iv) -mi iv) 



(C) -ni (C) h\'> (*x) 
-0-:(C) -C";(C) ^xh^U^x) 



0, (34) 



from which equations similar to Eq. pip are obtained. 
The roots of these equations gives the eigenvalues Eni in 
the range Ao < Eni < A2 that are obtained numerically. 

It is apparent from Eqs. ([29| and (p4| that i? = Ao is 
always the degenerate eigenvalues for / > 1 when A2 > Aq 
as ao = and ji (^ = 0) = for alH > 1 when £' = Ao. 
Also En,i=i = is the root of Eqs. ([29l) and ([Ml) as 
V''n,i=i — d(j)/dr has already satisfied the continuity con- 
dition Eq. (ITS)) that gives a part of the roots of Eqs. ([^^ 
and (1341). 



2. Spinodal regime 

In this case it is a simple textbook problem of a particle 
confined within a square well potential 



V (r) 



-|Ai|, 0<r<rB, 
A2, TB < r. 



(35) 



The radial solution for the bound state is now given by 
^ji (air) , < r < tb. 



Rni (r) 



Bh[^-' (ia2r) , tb < 



(36) 



The constants A and B will be determined from an equa- 
tion similar to Eqs. ((29)) and p4l) that leads to 



- ^Xhl'^' i^x) 31 (C) + CJ; (C) (ix) = 0. (37) 
This equation can be transformed into 



,(1) 



2c 



B^ 



(38) 



for I = 0^, which can be solved graphically. 

Since we are interested in the negative eigenvalues E < 
0, we consider the case when £' = in Eq. (|55)). which 
leads to 



(39) 
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as C = \/|Ai|/2crB and x = ^J^iI'^ctb when _B = 0. 
The roots of Eq. (pQ)) Co , Ci > ■ • ■ : gives the reduced ra- 
dius ^|Ai| /2crB for which the eigenvalues becomes zero 
{En.o — 0)i and are the lower bound for the appearance 
of one, two, • • • negative eigenvalues with I = Q. 
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three /=0 negative eigenvalues 


>^ — 'siii 




two /=0 negative eigenvalues 


^ 37U/2 


^, 


one /=0 negative eigenvalue 


■ jt 




^={\XM2cy'\ fromEq. (21) 


^"^ 7l/2~" 
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FIG. 4: The values Cn = /2crB ,n = 0, 1, 2, . . . which 

give the eigenvalues i?n,i=o = with I = 0. This eigenvalue 
En^i^o becomes negative if the matching radius vb satisfies 
(^n < •\/|Ai|/2crs. Also shown is the reduced matching radius 
Cb ~ \/\Xi\/2crB determined from the matching condition 
Eq. (|2ip of the stationary profile of nucleus. 



Figure m shows the roots Cn = y^\Xi\/2crB of Eq. (p9| . 
which gives the zero eigenvalues En^i=o — for I = 0, 
as the function of the ratio A2/IA1I. When A2 = 0, they 
are given by Co = 7r/2, Ci = 3^/2, C2 = 57r/2, . . 
Then, one negative eigenvalue with I = exist when 
the matching radius vb satisfies Ci > \/\^i\/2crB > Co, 
and two eigenvalues exist when C2 > ^|Ai|/2crs > Ci 
etc. As the barrier A2 makes the confinement more effec- 
tive, the number Cn increases as the function of A2/IA1I. 
Since Co is always smaller than the reduced matching ra- 
dius Cs — \/\M\/2crB determined from the matching 
condition Eq. (|2ip . we always have one negative eigen- 
value and, therefore, the lowest eigenvalues En=o.i=Q is 
always negative, which corresponds to the isotropically 
growing/shrinking spherical nucleusl'^ with / = 0. Other 
eigenvalues with I = are always positive as Cn > Cb for 
n=l,2,.... 

Similarly, Eq. (j37l) for I = 1 can be transformed into 

^-^, = ^ + 4 C' + X^-^-^rl, (40) 
C X X 2c 

Again, the condition for the zero eigenvalues E — with 
Z = 1 is given by 

It is easy to show that Eq. ([2T|) reduces to Eq. (|4T|). 
Therefore the reduced matching radius Cs = -^/jAi | /2crB 
determined from the matching condition Eq. (I^TI) of the 



stationary profile with the lowest free-energy always sat- 
isfies Eq. (HJ). Then, the eigenvalues for n — and 
Z = 1 is always zero (i?„=o,/=i — 0) since Co — Cb- 
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three /=1 negative eigenvalues 


in 




two /=1 negative eigenvalues 






one /=! negative eigenvalue 


f 7C 




= ilXMIc-yX fromEq. (21) 



0.5 1 1.5 2 

>.,/IX,l 

FIG. 5: The values Cn = ^|Ai| /2crB,n = 0,1,2,... which 
gives the zero eigenvalue E^.i^i = for 1 = 1, which becomes 
negative if Cn < \/\X\ \ /2crB- The reduced matching radius 
Cs = \/\Xi I /2crs determined from the matching condition 
Eq. (|2ip of the density profile of nucleus coincides with Co 
(Co = Cs)- 



Figure [5] shows the roots Cn — v^|Ai|/2crs of Eq. (|1T|) 
which gives the zero eigenvalue Enj=i ~ for ? = 1 as 
the function of the ratio A2/IA1I. When A2 — 0, they are 
given by tt, 2tt, Stt, . . .— . Again, the roots Cn increases 
as the function of A2/IA1I. Since the reduced matching 
radius Cs = \/\M\/2crB determined from the matching 
condition Eq. (PT|) coincides with Coj the lowest energy 
En=o 1=1 with n ~ and ^ = 1 is always zero {En=o 1=1 — 
0). ^ 

If the matching radius (Cb) is erroneously cho- 
sen from the multiple roots of Eq. (I?T|) such that the 
excited state with n = 1 and I = 1 has zero eigen- 
value En=i.i=i = (Ci = Cb), then not only the I = 
ground-state eigenvalues £'„=o,;=o but also the excited- 
state eigenvalues En=ij=o and £'„=o,i=i become nega- 
tive from Figs. |4] and [5] Therefore any numerical error 
in the determination of the stationary profile (j)s (r) of 
the critical nucleus could result in the multiple negative 
eigenvalues appearing. 

D. Numerical example and discussions 

In order to study the stability problem of the criti- 
cal bubble, we have solved Eqs. dM]) and to 
calculate the full spectrum of the bound-state eigenval- 
ues with / < 2 as the function of the scaled under- 
saturation Afj,/ Afispin- Figure [5] shows the bound-state 
eigenvalues E for the case (i) of Tab. ID Since Aq > A2 
(Tab. m, the square- well potential has the shape shown 
in Fig. Ela). The continuum states have the energy 
i?/|Ai| > A2/IA1I = 0.6. These continuum states de- 
scribe the capillary-like waves induced around the spher- 
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ical surface of nucleus. Since this continuum state starts 
from E = X2, there always exist minimum energy A2 to 
excite capillary wave. 



spherical growing mode as there will be an extra negative 
eigenvalue En=o,i=i with ^ = 1. 
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FIG. 6: Three bound-state eigenvalues with I — 0,1,2 as 
functions of the scaled under-saturation A/x/A/ispin for the 
case (i) of Tab. [H Only the ground state with / — has 
a negative eigenvalues. The lowest eigenvalues with i = 1 is 
always zero. There is only one negative eigenvalues with I = 
from the CNT regime to spinodal regime up to the spinodal 
point. 



There is only one bound state for each Z = 0,1,2 in 
Fig-El These eigenvalues depend on the under-saturation 
Afi in the CNT regime as the radius ta and tb of the 
potential well depend on the under-saturation (Fig. [2]) ■ 
Only the ground state eigenvalue for I — becomes neg- 
ative and decreases while the ground state eigenvalues 
for Z > 2 increase as the absolute magnitude |A/i| is in- 
creased toward the spinodal. A similar behavior of the 
eigenvalues with different I was observed in the numerical 
results of Varea and Robledoi. However, the eigenvalues 
become constant in the spinodal regime in our model as 
the matching radius tb becomes constant in this regime. 
Only the ground state with I — has the negative eigen- 
values that describes the growing/shrinking nucleuaii^ 
that preserves the spherical symmetry. This single neg- 
ative eigenvalue survives up to the spinodal point. The 
lowest eigenvalues with Z = 1 is always zero which means 
that the deformation with the form 



1 (r) oc 



cos 6 {d(j)s/dr) 
sin 6 {d(j)s/dr) 



(42) 



does not cost energy^ as Y;=i,m(6', ^) cos 6*, sin 6''^^. 
This deformation, in fact, describes the translation of 
the center of spherical nucleus'^^, and, therefore, deos not 
cost energy. 

It is possible to choose larger tb {(b) from Eq. (pij) 
that corresponds, for example, to Ci in Fig. [Sj How- 
ever, this false stationary profile has an artificial density 
oscillation that can be anticipated from Eq. (IT^ . and 
has a higher free-energy. Figures [5] and [S] also indicate 
that this false stationary state is unstable against non- 
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FIG. 7: The same as Fig.[6]for the case (ii) of Tab.|ll Now the 
bound state for I = 2 disappears before the spinodal regime 
is reached as the barrier A2 is low. 

Figure [7] shows the bound state energy for the case (ii) 
of Tab. m This is also the case shown in Fig. Ela). Now 
the continuum state starts from i?/|Ai| > 0.45. The I = 2 
bound state increases its energy as the under-saturation 
Afi is increased and finally it disappears into the contin- 
uum state as the confinement A2 is weak. Again there 
is only one negative eigenvalue with I = and one zero 
eigenvalue with I — 1 which can survive up to the spin- 
odal point. 

Figure [5] shows the bound state energy for the case (iii) 
of Table m In this case, the square- well potential has the 
shape shown in Fig. ^h) as A2 > Aq (Tab. |T|. In this 
case we can expect more complex energy diagram for the 
bound state. In particular, the degenerate bound states 
with E = Xq appear for all / > 1 at E/ |Ai| = 0.5 in 
the CNT regime. These bound states cannot survive in 
the spinodal regime as the inner barrier with E = Xq 
disappears in this regime. Once again there is only one 
negative eigenvalue with I = and one zero eigenvalue 
with Z = 1 up to the spinodal point. 

Therefore, the stationary profile of the critical bub- 
ble, which satisfies the matching condition of the den- 
sity profile Eqs. and (^(7)) and has the lowest free- 
energy, has only one unstable growing mode of the fluc- 
tuation with spherical symmetry. The time scale t of 
the growth or shrinkage is given by the inverse of the 
absolute magnitude \En=o,i=o\ of the negative eigenvalue 
with n = and / = from Eq. ITTI) . that is roughly 
given by r ~ 1/ (ri5„=o,z=o) 2/ (r|Ai|) in the spinodal 
regime (see Fig. [6] to [8|) . Therefore the curvature |Ai| at 
the top of the barrier of the free-energy Aut {(f) in Fig. [1] 
plays a crucial role in determining the timescale of how 
fast the saddle point is crossed once the critical bubble 
is formed. Of course, the time scale of nucleation that is 
the time necessary to form a critical bubble is given by 
the nculeation rate, which is determined by the work of 
formation of critical bubble^"— li^ii^ii^. 
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FIG. 8: The same as Fig.[6]for the case (iu) of Tab. HI In this 
case new bound-states with the energy E/ |Ai| = 0.5 appears 
for alU > 1 in the CNT regime. 

This conclusion persists up to the spinodal point. A 
diverging compressibility A2 ^ as Afj, A/igpin will 
not affect the timescale of nucleation even near the spin- 
odal. Also, any fluctuation other than the spherically 
growing mode that correspond to the negative eigenval- 
ues i?„=o,;=o with n = and I — stays stable up to the 
spinodal point. Therefore any structural anomaly22i2^iM 
near the spinodal will start during the growing stage af- 
ter the nucleation rather than at the nucleation stage of 
nucleus. Also this anomaly would be the effect of the in- 
teraction of multiple nuclei, such as coalescence or coars- 
ening^"— during the growing stage rather than the effect 
of the instability at the nucleation stage. Since we look at 
the stationary critical nucleus, the coupling of the order 
parameter to the diffusion and the heat flow is beyond 
the scope of the present work. These effects will also play 
crucial role in the structural anomaly of nucleus during 
the growing stage after the nucleatio n'^^'^^i'^^ 

IV. CONCLUSION 

In this study, a square-gradient density-functional 
model with a triple-parabolic free energy was used 
to study the stability of the critical bubble of ho- 
mogeneous bubble nucleation. By using this square- 
gradient model2.!i2 instead of the original density func- 
tional model^"— , the stability problem was reduced from 
the eigenvalue problem of the matrix with roughly 10'^ to 



10^ elements to the text-book problem of the eigenval- 
ues of the Schrodinger equation that describes a particle 
confined within a square-well potential. The negative 
eigenvalues of the bound state of the Schrodinger equa- 
tion represents the growing mode of the fluctuation and 
its magnitude determines the time scale of this grow- 
ing fluctuation. We found that there is only one nega- 
tive eigenvalue that corresponds to the spherically grow- 
ing/shrinking mode of nucleation. Our result confirmed 
the conclusion theoretically predicted^ and numerically 
obtained^ by other authors that the critical nucleus is lo- 
cated at the saddle point of the free-energy landscapei^. 

Furthermore, we have confirmed that this single nega- 
tive eigenvalue persists up to the spinodal point. There- 
fore no fractal or ramified structure^ is expected at the 
nucleation stage. Also, it is clear from our analysis that 
the stability analysis using the Schrodinger equation or 
the stability matrix is limited only for the critical nucleus 
at the saddle point in the free-energy landscape. In order 
to study the embryonic bubble before and after crossing 
the saddle point as the critical bubble, some authors^"— 
studied a constrained density functional model that were 
defined not with a simple smooth density profile but 
rather by a fixed particle number contained within a 
given small volume. It is not clear if the appearance of 
the negative eigenvalues for the stability matrix of this 
constrained system means the unstable growing mode of 
the growing bubble after crossing the saddle point. 

Finally, our stability analysis using the Schrodinger 
equation can be applicable to the special case of the den- 
sity functional model for the Yukawa fluid as the density 
functional for the Yukawa fluid can be transformed ex- 
actly into the square-gradient functional^i^. The sta- 
bility analysis of critical nucleus of this Yukawa model 
fluid will be presented elsewhere. 
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